Method and devices for image reconstruction

ABSTRACT

A method for reconstructing image data from x-ray data measured with an imaging system having at least one photon-counting detector includes obtaining a representation of data measured by the photon-counting detector. The method also includes generating first image data based on a projection based first functional using a first algorithm, the projection based first functional being dependent on the representation of data. The method also includes updating, based on a second functional that includes a model of at least one physical effect not included in the projection based first functional, the first image data to obtain second image data. The invention also provides an image processing device configured to reconstruct image data from x-ray data measured with an imaging system including at least one photon-counting detector as well as a corresponding computer program product.

TECHNICAL FIELD

The present disclosure relates to methods, devices and computer programsfor reconstructing image data based on performed x-ray measurements.

BACKGROUND

Reconstructing image data from photon-counting measurements may forexample include computing or estimating a vector a=(a_(jn)) of basiscoefficients of basis j=1, . . . , M_(e) and image pixel n=1, . . . ,M_(p). based on a vector of photon counts m=(m_(ik)) in energy bin i=1,. . . , M_(e) and detector element k=1, . . . , M_(d).

The photon counts may for example be obtained by direct readout pf aphoton-counting x-ray detector. In another example, the photon countsmay be the output of a post processing scheme operating on the countsobtained by readout from a photon-counting detector. The post-processingscheme may for example involve summation, filtering, averaging, andapplication of correction factors or correction terms.

In another example, the image data may consist of a single image,consisting of for example one of the basis coefficients or a combinationthereof. In yet another example, the image data may consists ofcoefficients in a non-pixelized representation of an image. For example,the elements of a may be Fourier coefficients, or wavelet coefficients,or coefficients in a representation of the image as a sum of blobs.

The vector a is, in a typical case, obtained through the optimization ofa function Φ(a, m) of a and m. A function which is optimized to find awill be referred to as a functional.

Such a functional, which is a function of measured data in theprojection domain, i.e. image data relating to the transmission throughor photon count values after an object, will be referred to as aprojection-based functional. A projection-based functional may be usedfor image-based material decomposition or reconstruction. In anotherexample, a projection-based functional may be a function of lineintegrals of the linear attenuation coefficient or of basis coefficientsalong a projection line.

A functional which is instead a function of measured data or datareconstructed from measured data, in the image domain, is called animage-based functional. Such a functional can be used to performimage-based material decomposition or reconstruction.

A functional Φ(a, m) may incorporate prior information about the imagedobject. For example, this prior information may be provided in the formof an edge-preserving regularizer, which penalizes rapid variations inthe image. In another example, the prior information may be provided inthe form of a discrepancy term penalizing a difference between thereconstructed image and a prior image.

As an example, in the absence of pulse pileup, i.e. at low photonfluence rates, the number of registered photons in energy bin k anddetector element k may be modelled as Poisson distributed with mean

$\begin{matrix}{\lambda_{ik} = {{\int\limits_{0}^{kVp}{{N_{0}(E)}{S_{i}(E)}{\exp\left( {- {\sum\limits_{j = 1}^{M_{b}}\;{{A_{jk}\left( {x,y} \right)}{\mu_{j}(E)}}}} \right)}{dE}}} + s_{ik}}} & (1)\end{matrix}$with A_(jk)=∫_(l) _(k) a_(j)dl integrated along projection ray k. Here,N₀(E) is the incident spectrum, S_(i)(E) is a weight function modelingthe sensitivity of the energy bin to different incident energy levelsand s_(ik) is the number of expected scatter counts from the imagedobject.

Eq. (1) does not model cross-talk between detector pixels, e.g. due tocharge sharing, K-fluorescence or internal Compton scatter, but this canbe included in the model by applying a linear operator B to the vector λof expected counts: λ^(c)=B λ. For example, B may be represented by asparse matrix with elements that are nonzero only fornearest-neighboring pixels. In another example, the elements of B maydecrease with increasing distance between pixels. In yet anotherexample, B may be obtained from a Monte Carlo simulation of a pencilbeam of photons impinging on one detector element and scattered intoneighboring detector elements.

Furthermore, Eq. (1) does not include pile-up. The effect of pile-up canfor example be modelled by letting the registered counts in oneprojection line be a function of the incident photons in all differentenergy bins in the same projection line: λ_(ik) ^(p)=f (λ_(1k) ^(c), . .. , λ_(M) _(e) _(,k) ^(c)) where λ_(ik) ^(p) is the number of registeredcount after pile-up in energy bin i and detector element k. For example,λ_(ik) ^(p) may be given by a paralyzable model or a non-paralyzabledetector model.

Modern CT reconstruction algorithms typically generate the reconstructedimage as a maximum a posteriori (MAP) estimate of the image given themeasured data. The MAP estimate may build on a complete model of therelationship between the registered counts and the image values, or itmay build on a simplified model of the relationship in order to simplifythe optimization algorithm. Compared to the analytical reconstructionalgorithms used previously, MAP reconstruction reduces noise and allowscorrection for detrimental effects such as scatter and optical blur. Forenergy-resolving, photon-counting CT, the MAP estimator of the vector ofbasis coefficients in the image a is, in a typical case:

$\begin{matrix}{\hat{a} = {{\underset{a}{argmin}{\sum\limits_{i = 1}^{M_{e}}\;{\sum\limits_{k = 1}^{M_{d}}\;\left( {{\lambda_{ik}^{p}(a)} - {m_{ik}{\lambda_{ik}^{p}(a)}}} \right)}}} + {R(a)}}} & (2)\end{matrix}$where λ_(ik) ^(p) and m_(ik) are the number of expected and registeredcounts, respectively, in energy bin i=1, . . . , M_(e) and detectorelement k=1, . . . , M_(d). λ_(ik) ^(p) includes the effect ofcross-talk and pileup is obtained from (1) with the blur operator B andthe pileup function ƒ applied. R(a) is an edge-preserving regularizer,which penalizes differences between neighboring detector elements. Theexpression that is optimized in (2) will be referred to as an MAPfunctional.

If the regularization term is not included, (2) becomes a maximumlikelihood (ML) estimator, and the functional to be minimized is calleda maximum likelihood functional

Since (2) is difficult to solve fast enough, it is common practice toreplace it with a simplified penalized weighted least squares estimator.This estimator is given by minimization of a penalized weighted leastsquares functional, for example:

$\begin{matrix}{\hat{a} = {{\underset{a}{argmin}{\sum\limits_{j = 1}^{M_{b}}\;{\sum\limits_{k = 1}^{M_{d}}\;\frac{\left( {{Ta} - {\hat{A}}_{jk}} \right)^{2}}{\sigma^{2}\left( {\hat{A}}_{jk} \right)}}}} + {R(a)}}} & (3)\end{matrix}$

Here, T denotes the forward ray transform operator and Â_(jk) is anestimate of the line integral A_(jk)=∫_(l) _(k) a_(j)dl along projectionray k. σ²(Â_(jk)) is the variance of Â_(jk). Â_(jk) can be obtained fromthe measured counts m_(ik) for each individual detector element usingmaximum likelihood estimation or a look-up table. Eq. (3) can becomputed quickly using e.g. the iterative coordinate descent (ICD)method or the separable quadratic surrogates (SQS) method, but it givesinferior image quality compared to (2) since it builds on a simplifiednoise model and ignores detector cross-talk and object scatter. Inparticular, (3) is based on modelling the noise as Gaussian instead ofPoisson, which is an approximation. There is thus a need for analgorithm that combines good image quality with fast computation time.

The publication “Multi-Material Decomposition Using Statistical ImageReconstruction for Spectral CT” by Y. Long and J. Fessler, IEEE Trans.Med. Imag. 33, pp. 1614-1626 (2014), relates to an iterativereconstruction method for penalized-likelihood reconstruction ofmaterial basis images. This method is initialized by a set of basisimages obtained through an image-domain material decomposition method.

U.S. Pat. No. 5,390,258 relates to a method of acquiring an image froman object, wherein a set of training images is used to generate aconvergent series expansion, and wherein the measured signals are usedto generate a truncated series expansion of an image of the object.

U.S. Pat. No. 6,754,298 relates to a method based on a statistical modelfor reconstructing images from transmission measurements with highenergy diversity.

U.S. Pat. No. 7,551,708 relates to a method of reconstructing materialdecomposed images from data from energy discriminating computedtomography detectors using the iterative coordinate descent (ICD)algorithm.

U.S. Pat. No. 9,165,384 relates to a method of image reconstructionwhich reconstructs a plurality of final component images of an objectbased on spectral projection data, wherein intermediate images are usedin the reconstruction algorithm and wherein correlations between theseintermediate images are taken into account in the algorithm.

U.S. Pat. No. 8,923,583 relates to a tomographic reconstruction methodwherein material component images are reconstructed by optimizing ajoint likelihood functional, which includes information on correlationsbetween the component sinograms.

US patent application 2016120493A1 relates to an x-ray CT imageprocessing method wherein a joint posterior distribution, based on aprior probability distribution, is used to estimate an x-ray absorptioncoefficient from measurements with different wavelengths.

U.S. Pat. No. 8,929,508 relates to a method of computing line integralsof basis coefficients through an object from x-ray photon transmissionmeasurements by computing a first approximation to the line integralsand combining the first approximation with a correction computed from acalibration phantom.

U.S. Pat. No. 6,907,102 relates to a method of image iterativereconstruction wherein a cross-section reconstruction vectorapproximately matching the projection data is determined using acomputed tomography model.

U.S. Pat. No. 7,885,371 relates to a method of tomographic imagereconstruction wherein a first reconstruction method, which convergesfaster on low spatial frequencies than on high spatial frequencies, isfollowed by a second reconstruction method which converges faster onhigh spatial frequencies than on low spatial frequencies.

U.S. Pat. No. 9,508,163 relates to a method of iterative tomographicreconstruction wherein each iteration of an outer loop includesiterative processing of an inner loop.

U.S. Pat. No. 9,020,230 relates to a reconstruction method employing anouter iteration loop and an inner iteration loop, wherein the inneriteration loop calculates a preconditioner used by the outer loop.

U.S. Pat. No. 6,256,367 relates to a method of correcting for artifactsdue to scatter CT images by Monte Carlo simulating photon scatter andsubtracting the simulated photon energy from the measured projectiondata.

There is, however, an ongoing need for a fast reconstruction algorithmthat gives good image quality.

SUMMARY

It is an object of the proposed technology to provide methods andcorresponding devices and computer programs that yields reconstructedimages of high quality fast.

It is a particular object to provide a method for reconstructing imagedata based on measurements performed by an x-ray system.

It is another particular object to provide an image processing devicefor reconstructing image data based on measurements performed by anx-ray system.

It is yet another object to provide a computer program forreconstructing image data based on measurements performed by an x-raysystem.

According to a first aspect there is provided a method of reconstructingimage data from x-ray data measured with an imaging system comprising atleast one photon-counting detector. The method comprises the step ofobtaining a representation of data measured by the photon-countingdetector. The method also comprises the step of generating first imagedata based on a projection based first functional using a firstalgorithm, the projection based first functional being dependent on therepresentation of data. The method further comprises the step ofupdating, based on a second functional that includes a model of at leastone physical effect not included in the projection based firstfunctional, the first image data to obtain second image data.

According to a second aspect there is provided an image processingdevice for reconstructing image data from x-ray data measured with animaging system comprising at least one photon-counting detector. Theimage processing device is configured to obtain a representation of datameasured by the photon-counting detector. The image processing device isalso configured to generate first image data based on a projection basedfirst functional using a first algorithm, the projection based firstfunctional being dependent on the representation of data. The imageprocessing device is also configured to update, based on a secondfunctional that includes a model of at least one physical effect notincluded in the projection based first functional, the first image datato obtain second image data.

According to a third aspect there is provided a computer program adaptedto reconstruct image data from x-ray data measured with an imagingsystem comprising at least one photon-counting detector and comprisesinstructions, which when executed by at least one processor, cause theprocessor(s) to:

-   -   read a representation of data measured by a photon-counting        detector of an x-ray imaging system;    -   generate first image data based on a first functional using a        first algorithm, the first functional being dependent on the        representation of data;    -   update, based on a second functional that includes a model of at        least one physical effect not included in the first functional,        the first image data to obtain second image data.

According to a fourth aspect there is provided a computer programproduct comprising the computer program of the fourth aspect.

The proposed technology provides for a fast algorithm that yields highquality image data. The provided algorithm has the virtue of being ableto generate an image that takes several physical effects into account,such as non-Gaussian statistics, detector cross-talk, pile-up andoptical blur. The provided algorithm is also fast since it builds onsolving a simplified optimization problem as a first step and thencorrecting the resulting first image with a few computationallyinexpensive steps.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a schematic diagram illustrating an x-ray imaging systemaccording to prior art.

FIG. 2 is a schematic diagram illustrating an alternative version of anx-ray imaging system according to prior art.

FIG. 3 is a schematic flow diagram illustrating a method forreconstructing image data according to the proposed technology.

FIG. 4 is a schematic flow diagram illustrating a particular version ofthe proposed method.

FIG. 5 is a schematic flow diagram illustrating another particularversion of the proposed method.

FIG. 6 is a schematic flow diagram illustrating still another particularversion of the proposed method.

FIG. 7 is a diagram illustrating a computer program implementationaccording to the proposed technology.

FIG. 8 is a diagram illustrating a particular embodiment of an imageprocessing device according to the proposed technology.

DETAILED DESCRIPTION

It may be useful to begin with a brief overview of an illustrativeoverall x-ray imaging system, with reference to FIG. 1. In thisnon-limiting example, the x-ray imaging system 100 basically comprisesan x-ray source 10, an x-ray detector system 20 and an associated imageprocessing device 30. In general, the x-ray detector system 20 isconfigured for registering radiation from the x-ray source 10 that mayhave been focused by optional x-ray optics and passed an object orsubject or part thereof. The x-ray detector system 20 is connectable tothe image processing device 30 via suitable analog processing andread-out electronics (which may be integrated in the x-ray detectorsystem 20) to enable image processing and/or image reconstruction by theimage processing device 30.

As illustrated in FIG. 2, another example of an x-ray imaging system 100comprises an x-ray source 10, which emits x-rays; an x-ray detectorsystem 20, which detects the x-rays after they have passed through theobject; analog processing circuitry 25, which processes the rawelectrical signal from the detector and digitizes it; digital processingcircuitry 40 which may carry out further processing operations on themeasured data such as applying corrections, storing it temporarily, orfiltering; and a computer 50 which stores the processed data and mayperform further post-processing and/or image reconstruction.

A challenge for x-ray imaging detectors is to extract maximuminformation from the detected x-rays to provide input to an image of anobject or subject where the object or subject is depicted in terms ofdensity, composition and structure. It is still common to usefilm-screen as detector but most commonly the detectors today provide adigital image.

Modern x-ray detectors normally need to convert the incident x-rays intoelectrons, this typically takes place through photo absorption orthrough Compton interaction and the resulting electrons are usuallycreating secondary visible light until its energy is lost and this lightis in turn detected by a photo-sensitive material. There are alsodetectors, which are based on semiconductors and in this case theelectrons created by the x-ray are creating electric charge in terms ofelectron-hole pairs which are collected through an applied electricfield.

In general, the x-ray photons, including also photons after Comptonscattering, are converted to electron-hole pairs inside a semiconductordetector, where the number of electron-hole pairs is generallyproportional to the photon energy. The electrons and holes are thendrifting towards the detector electrodes, then leaving the detector.During this drift, the electrons and holes induce an electrical currentin the electrode, a current which may be measured, e.g. through a ChargeSensitive Amplifier (CSA), followed by a Shaping Filter (SF).

As the number of electrons and holes from one x-ray event isproportional to the x-ray energy, the total charge in one inducedcurrent pulse is proportional to this energy. The current pulse isamplified in the CSA and then filtered by the SF filter. By choosing anappropriate shaping time of the SF filter, the pulse amplitude afterfiltering is proportional to the total charge in the current pulse, andtherefore proportional to the x-ray energy. Following the SF filter, thepulse amplitude is measured by comparing its value with one or severalthreshold values (Thr) in one or more comparators (COMP), and countersare introduced by which the number of cases when a pulse is larger thanthe threshold value may be recorded. In this way it is possible to countand/or record the number of X-ray photons with an energy exceeding anenergy corresponding to respective threshold value (Thr) which has beendetected within a certain time frame.

When using several different threshold values, a so-calledenergy-discriminating detector is obtained, in which the detectedphotons can be sorted into energy bins corresponding to the variousthreshold values. Sometimes, this type of detector is also referred toas a multi-bin detector.

In general, the energy information allows for new kinds of images to becreated, where new information is available and image artifacts inherentto conventional technology can be removed.

In other words, for an energy-discriminating detector, the pulse heightsare compared to a number of programmable thresholds in the comparatorsand classified according to pulse-height, which in turn is proportionalto energy.

Having briefly described how photon counting detectors work and how themeasurements may be used to reconstruct images, in what follows therewill be described methods and devices that provides improved imagereconstruction features.

The inventors have appreciated that the optimum of a first, simplifiedfunctional, for example a penalized weighted least squares functional,may be a good approximation to the optimum of a second, more complicatedfunctional.

The inventors have further appreciated that a projection-basedfunctional is able to model the image acquisition more accurately thanan image-based functional. Therefore, it is preferable to let the first,simplified functional, be a penalized weighted least squares functional.

The publication “Multi-Material Decomposition Using Statistical ImageReconstruction for Spectral CT” by Y. Long and J. Fessler, IEEE Trans.Med. Imag. 33, pp. 1614-1626 (2014), relates to a method of using a setof basis images resulting from image-based material decomposition toinitialize an optimization of a functional which contains a morecomplete model of the imaging system, including a model of Poissonnoise. However, the method of the publication does not include usingbasis images resulting from the optimization of a projection-basedfunctional to initialize the optimization.

We disclose a novel method of solving the image reconstruction problem.To this end there is provided a method of reconstructing image data fromx-ray data measured with an imaging system comprising at least onephoton-counting detector. The method comprises the step of obtaining S1a representation of data measured by the photon-counting detector. Themethod also comprises the step of generating S2 first image data basedon a projection based first functional using a first algorithm, theprojection based first functional being dependent on the representationof data. The method further comprises the step S3 of updating, based ona second functional that includes a model of at least one physicaleffect not included in the projection based first functional, the firstimage data to obtain second image data. The method is schematicallyillustrated in the flow diagram of FIG. 3.

The proposed method may for example be based on first solving a simplerproblem, e.g., the problem of solving the functional defined by formula(3) above, by using a first optimization algorithm and then applying atleast one update so that the solution becomes a better approximation ofthe full MAP problem defined by formula (2) above. In general, it buildson obtaining first image data by optimizing a first, simplifiedfunctional and subsequently applying one or more update steps so thatthe resulting image data is a better approximation than the first imagedata to the optimum of a second functional, which models additionalphysical effects not modelled by the first functional.

Put in slightly different words, the proposed technology provides amethod for reconstructing images based on measurements performed by aphoton counting detector. The measurements are provided, in a step S1 asinput in order to be able to generate, in a step S2, image data based ona first projection based functional. The measurement data is provided ina representation that is suitable for computing the first functional.The measurement data may thus be provided in a pre-processed form, i.e.the raw measurement data may have been subject to pre-processing orpre-handling whereby a particular measurement data representation isobtained. The particular representation should be chosen so as to obtaina form that is suitable for use as input in the first functional inorder to generate the first image data. The measurement data may inparticular be tomographic x-ray data measured with an imaging systemcomprising at least one photon-counting detector. The method alsocomprises to update the generated first image data in order to obtainsecond image data. The update may be done in several ways which will bedescribed below. The updates are however all done based on a secondfunctional that provides a more detailed model of the physical effectspresent during x-ray detection. A particular example illustrating thedifferences between the first functional and a second, more detailed,functional is given by the functional defined by formula (3) and thefunctional defined by formula (2) provided earlier.

According to a particular embodiment of the proposed method, the step S2of generating first image data might comprise using a first algorithm onthe first functional where the first functional is a penalized weightedleast squares functional.

According to another particular embodiment of the proposed method, thestep S2 of generating first image data comprises using a first algorithmon the first functional where the first functional is a maximumlikelihood or a maximum a posteriori functional.

Still another embodiment provides a method wherein the first algorithm,also referred to as a first optimization algorithm, is a Newton method,an iterative coordinate descent method, a separable quadratic surrogatesmethod, an expectation maximization method, a conjugate gradient method,an ordered subset method, or a combination of a plurality of themethods.

In this way, the optimum of the simplified functional, which is easy tocompute, can be corrected in one or a small number of update steps, toyield an approximate optimum of the second functional. In this way, alarge number of iterations of an optimization algorithm converging tothe optimum of the second functional, which could be computationallyexpensive, is avoided.

Yet another embodiment provides a method where the update is performedin a first manner utilizing a particular term from the secondfunctional. This embodiment provides a method wherein the step S3 ofupdating the first image data to obtain second image data comprises toperform an image updating algorithm on the first functional when atleast one term from the second functional has been added to the firstfunctional, the at least one term providing a model of at least onephysical effect not included in the first functional.

Stated slightly differently, the first functional may be used in itsoriginal form, see e.g. the functional representation provided byformula (3) above, to yield first image data a. This first image data amay then be updated by adding a particular term from the secondfunctional, see e.g. functional representation provided by formula (3)for a particular second functional, that provides a model of a physicaleffect not included in the first functional. This particular term may beseen as a perturbation of the first functional, i.e. the firstfunctional is perturbed by the particular term from the secondfunctional. In this way it will be possible to update the first imagedata a by inserting it into the perturbed first functional and performan image updating algorithm in order to yield second image data a*. Thephysical effects modelled by the particular term may include a model ofthe Poisson noise, pile-up, optical blur due to the finite size of thedetector and focal spot, or detector cross-talk. It is possible to addseveral distinct terms from the second functional to the firstfunctional. The different terms may model several different physicaleffects. The image updating algorithm used for the initial firstfunctional, e.g. the non-perturbed first functional, may be a differentimage updating algorithm to the image updating algorithm that is used toobtain updated image data by applying it to the perturbed firstfunctional. This embodiment provides a way for updating the image databy using a particular term, or particular terms, from the secondfunctional.

The disclosed method works best if the second functional can be regardedas a small perturbation of the simplified functional. This is often thecase if the second functional, for example, is based on a representationof one or several physical effects that are not modelled by thesimplified functional. The physical effects may include a model of thePoisson noise, pile-up, optical blur due to the finite size of thedetector and focal spot, or detector cross-talk.

According to a possible embodiment of the proposed technology there isprovided a method where said image updating algorithm is performedutilizing a series expansion of the optimum of the perturbed firstfunctional around the first image data with respect to at least oneparameter determining the magnitude of at least one particular term fromthe second functional. To this end there is provided a method whereinthe step of updating the first image data comprises to perform a seriesexpansion of at the optimum of the perturbed first functional around thefirst image data with respect to at least one parameter describing themagnitude of at least one particular term in the second functional andcompute, based on the series expansion, updated image data where atleast one term in the series expansion of the optimum has been added tothe first image data to thereby obtain second image data. Preferably theupdate comprises the addition of at least one term in a series expansionthat converges to an optimum of the second functional.

In different words, beside the fact that it is possible to add severaldistinct terms from the second functional to the first functional, it isalso possible to add one or several image data correction terms obtainedfrom a series expansion representing the influence of one or several ofthe terms from the second functional on the optimum of the perturbedfirst functional, to the first image data. Calculation of the image datacorrection terms from the series expansion may be done through directapplication of an analytical formula, or by solving an equation by aniterative method. The series expansion may for example be a Taylorexpansion around a perturbation strength parameter d.

In the present disclosure an optimum of a functional refers to either amaximum point or a minimum point. A maximum point of a functional L isan input data vector a of the functional L such that the value L(a) ofthe functional evaluated with the data vector as input is larger thanthe value L(a′) of the functional evaluated for any other input datavector a′ in a set of permissible input data vectors. Similarly, aminimum point of a functional L is an input data vector a of thefunctional L such that the value L(a) of the functional evaluated withthe data vector as input is smaller than the value L(a′) of thefunctional evaluated for any other input data vector a′ in a set ofpermissible input data vectors. There are two types of functionals usedin image reconstruction: functionals that should be maximized andfunctionals that should be minimized. Which of these two categories agiven functional belongs to can be seen on the mathematical formulationof the functional. The two types of functionals are closely related,since a functional that should be maximized can be turned into afunctional that should be minimized by multiplying it with −1, and viceversa. Optimizing a functional amounts to either maximizing it, in caseit is a functional that should be maximized, or minimizing it, in caseit is a functional that should be minimized. In practice, the optimummay not be possible to calculate exactly, so that optimization in apractical situation may amount to calculating an approximation to theminimum point or to the maximum point.

By way of example, if Φ(a, m) denotes the target functional to beminimized in (2) and Φ⁰ (a, m) the target functional in (3), Φ(a, m) canbe expressed as Φ⁰ (a, m)+dΦ¹(a, m) where d is a perturbation strengthparameter. If a=h(d, m) denotes the mapping from registered data m to anestimated set of basis images a, h(d, m) can be expanded in a Taylorseries in d around d=0. For example, to first order in d, h(d, m)=h⁰(m)+d·h¹ (m), where the image correction term h¹ is given byH _(Φ) ₀ h ¹(m)=∇_(a)Φ¹(a,m)|_(a=h) ₀ ₍ m)  (4)where H_(Φ) ₀ is the Hessian matrix of Φ⁰ with respect to a. Solvingthis quadratic equation for h¹ typically requires an amount ofcomputational power comparable to two iterations with an iterativemethod for solving (3). Since direct iterative optimization of (2) mayrequire hundreds of iterations, this can be very time-saving. Ifnecessary, higher order corrections can be calculated similarly.

It should be noted that the above mentioned perturbation strengthparameter d in certain embodiments may be set to one after the Taylorexpansion has been performed. In particular cases it is however notessential that the computation of the first functional, when perturbedby terms from the second functional, yields a converging solution. Thefirst few terms that are computed may yield a good approximation. In thelatter cases it may not be necessary to set d to one.

For example, the proposed method can combine the ability of the full MAPmethod (2) to model Poisson noise, detector cross-talk and objectscatter with the speed of the penalized weighted least squares problem(3). The fact that a well-established reconstruction method may be usedas a part of the proposed method also facilitates its introduction inclinical CT scanners.

The second algorithm referred to above may optionally be a Newtonmethod, an iterative coordinate descent method, a separable quadraticsurrogates method, an expectation maximization method, a conjugategradient method, an ordered subset method, or a combination of aplurality of the methods.

In another aspect of the invention, the image data can be updated inseveral steps, such that the result of a first set of updates with afirst update strategy is used as input to a second set of updates with asecond update strategy. The above proposed algorithms may thereby becombined sequentially.

According to a particular embodiment of the proposed technology there isprovided a method that can be used to select particular image data thatcan be used as reconstructed image data. To this end there is provided amethod that also comprises the steps of:

computing S4 the value of the second functional when the secondfunctional:

a) depends on the first image data, and

b) depends on the second image data; and

comparing S5 the computed values, and

selecting S6, based on the comparison, image data to be used asreconstructed image data.

One purpose of this particular embodiment is to select the particularimage data that best approximates the optimum of the second functional.A particular criteria that can be used where e.g. the second image datais selected if it is considered a better approximation to, e.g., theoptimum of the second functional than the first image data.

Consider the following non-limiting example of the proposed embodiment.At first a representation of the measured data is obtained in a step S1.In a step S2 first image data is generated based on a projection basedfirst functional, denoted L1 using a first algorithm and the obtainedrepresentation of data. At this point in the procedure first data,referred to as a, is obtained. The obtained data, a, may now be updatedbased on the second functional, referred to as L2. The second functionalmay be written asL2=L1+ΔL,Where the first functional L1 does not contain the terms used asperturbations for the first functional, i.e. the perturbation terms arecomprised in ΔL. ΔL includes terms that model of at least one physicaleffect that is not included in the projection based first functional. Asdescribed earlier, there are now several distinct ways to update thefirst image data in order to obtain second image data. One particularway is to add one or several terms comprised in ΔL to the firstfunctional in order to get L1+ΔL. Instead of including the completedifference between the first functional and the second functional in thecorrection term ΔL, one possibility is to expand one or several of theterms comprised in ΔL in a Taylor series and add at least one of theterms to L1.

This perturbed first functional may now be subject to an algorithm,either different from the one used to generate first image data or thesame algorithm, whereby second image data, referred to as a*, isobtained by updating the first image data in the step S3. This secondalgorithm may exemplarily be an iterative algorithm. In another aspectof the invention, the second algorithm may include Taylor expanding theoptimum of the perturbed first functional, viewed as a function of theinput data, around the first image data a. The Taylor expansion isthereby made with respect to a parameter d determining the magnitude ofthe perturbation term. From this Taylor expansion, one or severalcorrection term can be computed, which are applied to a and therebyyield a*.

Regardless of which particular way that is used to incorporate theperturbation/term, second image data will be obtained by applying analgorithm to the perturbed first functional. At this stage in theprocedure first and second image data, a and a*, is obtained. The methodnow proceeds with the aim of determining which particular image data ofthe first and second image data that best approximates the secondfunctional. That is, the method proceeds and determines which of thefirst and second image data that best approximates the complete solutionto the second functional. This may also be seen as determining which ofthe first and second image data that converges to a solution to thesecond functional. This may be obtained by first computing, in a stepS4, the value of the second functional L2 when the second functionaldepends on the first image data, i.e. L2 (a) and when it depends on thesecond image data, i.e. L2 (a*). After this the computed values L2 (a)and L2 (a*) are compared in a step S5 in order to select, in a step S6,image data to be used as reconstructed image data based on thecomparison. It should here be noted that a functional that is evaluatedbased on a particular image data, which might be a vector or even avector valued function, provides a scalar as output. The purpose of thecomparing step and selecting step is to select the particular image datathat best approximates a solution to the second functional. To clarifywhat is meant by a better approximation in the present case, a firstimage data vector a₁ is a better approximation than a second image datavector a₂ to a third image data vector a₃ if a mathematical measure ofthe discrepancy between a₁ and a₂ is less than the mathematical measureof the discrepancy between a₂ and a₃. A mathematical measure of thediscrepancy between a and a′ may for example be an Euclidian norm ofa−a′, or a weighted Euclidean norm of a−a′, or an L^(p) norm of a−a theabsolute value of the largest element of a−a′, or the Kullback-Leiblerdivergence between a−a′.

According to a particular version of the above mentioned embodiment theselected image data may be used as input in an updating scheme forupdating the reconstructed image data.

According to this aspect of the invention, the image data can be updatedin several steps, such that the result of a first set of updates with afirst update strategy is used as input to a second set of updates with asecond update strategy. The above proposed algorithms may thereby becombined sequentially.

An optional embodiment of the methods described above provides a methodwherein the step of computing the value of the second functionalcomprises to compute an estimate of the value second functional byperforming at least one step in an algorithm.

That is, the first image data and the second image data may be fed intothe second functional and be subjected to one or several steps in analgorithm in order to yield an approximation of the value of the secondfunctional. The particular image data to be used as reconstructed imagedata may then be selected, in step S6, based on the obtainedapproximation.

A particular version of the proposed embodiment provides a wherein thestep S6 of selecting the image data to be used as input comprises toselect the image data that yields the lowest or highest value for thesecond functional. With the value of the second functional is alsointended an approximate value of the second functional that may havebeen obtained by performing one or several steps of an algorithm appliedto the second functional.

In another embodiment of the invention, the at least one update of thefirst image data may be one or several iteration steps of a secondoptimization algorithm, such that the output of the algorithm convergesto the optimum of the second functional. For example, the steps of thesecond optimization algorithm may be more computationally heavy than thefirst algorithm, so that it is preferable to use the first algorithm tooptimize the first functional and one or a few steps of the secondalgorithm to update the image data, thereby approximating the optimum ofthe second functional.

For example, the second optimization algorithm may be a Newton method,which is may be computationally expensive but still possible to use fora small number of iterations.

In a preferred embodiment, the at least one update comprises the firstpart of a systematic update scheme the result of which, if the updatesare applied repeatedly, converges to the solution of the secondfunctional. With a systematic update scheme is meant a set of ruleswhich allows calculation of a series of updated sets of image data, froma first set of image data and, possibly, also from the result ofapplying previous updates to a first set of image data.

Such an update scheme is different from ad-hoc corrections since it willconverge to the statistically optimal image estimator provided that thesecond functional, for example the full problem (2) can be regarded as asmall perturbation of the first functional, for example (3).

The physical effect modelled by the second functional may in all of thedescribed embodiments relate to one of the following, or a combinationthereof: Poisson noise statistics, optical blur, pile-up, detectorelement cross-talk and object scatter.

According to a preferred embodiment of the proposed method there isprovided a method wherein the first image data comprises at least onebasis material image.

This particular embodiment utilizes a technique enabled byenergy-resolved x-ray imaging. That is, a technique commonly referred toas basis material decomposition. This technique utilizes the fact thatall substances built up from elements with low atomic number, such ashuman tissue, have linear attenuation coefficients μ(E) whose energydependence can be expressed, to a good approximation, as a linearcombination of two (or more) basis functions:μ(E)=a ₁ƒ₁(E)+a ₂ƒ₂(E).where ƒ_(i) are the basis functions and a_(i) are the correspondingbasis coefficients. If there is one or more element in the imaged volumewith high atomic number, high enough for a k-absorption edge to bepresent in the energy range used for the imaging, one basis functionmust be added for each such element. In the field of medical imaging,such k-edge elements can typically be iodine or gadolinium, substancesthat are used as contrast agents.

Basis material decomposition has been described in Alvarez and Macovski,“Energy-selective reconstructions in X-ray computerised tomography”,Phys. Med. Biol. 21, 733. In basis material decomposition, the lineintegral A_(i) of each of the basis coefficients a_(i) is inferred fromthe measured data in each projection ray l from the source to a detectorelement. The line integral A_(i) can be expressed as:A _(i)=∫_(l) a _(i) dl for i=1, . . . , N,where N is the number of basis functions. In one implementation, basismaterial decomposition is accomplished by first expressing the expectedregistered number of counts in each energy bin as a function of A_(i).Typically, such a function may take the form:

$\lambda_{i} = {\int\limits_{E = 0}^{\infty}{{S_{i}(E)}{\exp\left( {- {\sum\limits_{j = 1}^{N}\;{A_{j}{f_{j}(E)}}}} \right)}{dE}}}$Here, λ_(i) is the expected number of counts in energy bin i, E is theenergy, S_(i) is a response function which depends on the spectrum shapeincident on the imaged object, the quantum efficiency of the detectorand the sensitivity of energy bin i to x-rays with energy E. Even thoughthe term “energy bin” is most commonly used for photon-countingdetectors, this formula can also describe other energy resolving x-raysystems such as multi-layer detectors or kVp switching sources.

Then, the maximum likelihood method may be used to estimate A_(i) underthe assumption that the number of counts in each bin is a Poissondistributed random variable. This is accomplished by minimizing thenegative log-likelihood function, see Roessl and Proksa, K-edge imagingin x-ray computed tomography using multi-bin photon counting detectors,Phys. Med. Biol. 52 (2007), 4679-4696:

${\hat{A}}_{1},\ldots,{{\hat{A}}_{N} = {{\underset{A_{1},\ldots,A_{N}}{argmin}{\sum\limits_{i = 1}^{M_{b}}\;{\lambda_{i}\left( {A_{1},\ldots,A_{N}} \right)}}} - {m_{i}\ln\;{\lambda_{i}\left( {A_{1},\ldots,A_{N}} \right)}}}}$where m_(i) is the number of measured counts in energy bin i and M_(b)is the number of energy bins.

From the line integrals A, a tomographic reconstruction to obtain thebasis coefficients a_(i) may be performed. This procedural step may beregarded as a separate tomographic reconstruction, or may alternativelybe seen as part of the overall basis decomposition.

When the resulting estimated basis coefficient line integral Â_(i) foreach projection line is arranged into an image matrix, the result is amaterial specific projection image, also called a basis image, for eachbasis i. This basis image can either be viewed directly (e.g. inprojection x-ray imaging) or taken as input to a reconstructionalgorithm to form maps of basis coefficients a_(i) inside the object(e.g. in CT). Anyway, the result of a basis decomposition can beregarded as one or more basis image representations, such as the basiscoefficient line integrals or the basis coefficients themselves.

According to a particular example of the proposed technology the firstimage data may be obtained by applying a first optimization algorithm onthe first functional. The first optimization algorithm is a Newtonmethod, an iterative coordinate descent method, a separable quadraticsurrogates method, an expectation maximization method, a conjugategradient method, an ordered subset method, or a combination of aplurality of the methods. The update of the first image data to obtainsecond image data may in the same example be obtained by the applicationof at least one iteration of a second iterative method optimizing thesecond functional. The second iterative method is a Newton method, aniterative coordinate descent method, a separable quadratic surrogatesmethod, an expectation maximization method, a conjugate gradient method,an ordered subset method, or a combination of a plurality of themethods.

The update of the first image data may moreover comprise the addition ofat least one term in a series expansion converging to an optimum of thesecond functional.

In another example the second optimization method may be an iterativecoordinate descent method, a separable quadratic surrogates method (alsocalled a separable paraboloidal surrogates method), an expectationmaximization method, a conjugate gradient method, an ordered subsetmethod, or a combination of these.

In another embodiment of the invention, the full Poisson noise model canbe included in the optimization problem that is solved iteratively,using for example a pixel-wise separable quadratic surrogates (PWSQS)method similar to that described in Y. Long and J. Fessler,“Multi-Material Decomposition Using Statistical Image Reconstruction forSpectral CT”, IEEE Trans. Med. Imag. 33, pp. 1614-1626 (2014). This maybe beneficial in situations where the Poisson noise model differs somuch from the Gaussian noise model that an inconveniently large numberof updates to the first image data are needed. The series expansionmethod is then used to correct for, for example, detector cross-talk,pile-up and object scatter.

In yet another embodiment of the invention, it may be that the detectorcross-talk alters the image so much that a few terms in the seriesexpansion are insufficient to correct for it. In this case, the detectorcross-talk must be included in the optimization of the first functionalby including a blur operator in (3), while the series expansion is usedfor taking the Poisson noise model and object scatter into account. Theinitial optimization may for example be made using an optimizationalgorithm based on the separable quadratic surrogates method.

In another embodiment of the invention, both Poisson noise and detectorcross-talk may be represented in the optimization of the firstfunctional, whereas object scatter, which is too expensive to compute ineach iteration, can be corrected for using the series expansion method.

To facilitate the understanding of the propose technology below willfollow some illustrative examples of the described method.

In the first example reference is made to FIG. 4, where a flowchart ofan exemplary implementation of an iterative reconstruction algorithm isshown. In each iteration, the measured data, or a processed version ofthe measured data, is compared with the forward projection of thecurrent image estimate. Based on the discrepancy between these two setsof projection data, a back-projection step is used together with priorinformation about the expected content of the image to calculate animage update which is applied to the current image estimate to form anew image estimate. This procedure is iterated until a stoppingcriterion is fulfilled. Such an iterative reconstruction algorithm mayexemplarily be used in the present invention for generating S2 firstimage data or for updating S3 the first image data to obtain secondimage data.

In the second example reference is made to FIG. 5, where a flowchart ofan exemplary implementation of the present invention is shown. Measureddata is obtained from a photon-counting detector and processed to yieldpreprocessed data. Said preprocessed data is then used as input to afirst reconstruction algorithm which is based on a first functional andgenerates a first reconstructed image as output. The first reconstructedimage is then updated in an image update step, where the update is basedon a second functional, thereby yielding a final reconstructed image.

In the third example reference is made to FIG. 6, where a flowchart of afurther exemplary implementation of the present invention is shown.Measured data is obtained from a photon-counting detector and processedto yield preprocessed data. The preprocessed data is then used as inputto a first reconstruction algorithm which optimizes a first functionalusing an iterative algorithm and generates first reconstructed imagedata as output. Based on the first reconstructed image data and based ona second functional, which includes a model of one or more physicaleffects not included in the first functional, one or more imagecorrections is calculated. The image corrections are then applied to thefirst reconstructed image data to yield final reconstructed image data.

Having described various embodiments of the proposed method in detail,in what follows we will describe an image processing device configuredto perform the proposed method. All advantages associated with theearlier describe method are applicable also to the image processingdevice.

Reference is now made to FIG. 1, where it is illustrated how an x-rayimaging system 100 comprising an image processing device 30 connected toan x-ray system 20 adapted for detecting x-rays transmitted from anx-ray source 20. The x-ray system provides the image processing device30 with representations of the measured data. The representations may beany representation that is suitable as input in an image reconstructionprocess. The particular ways and representations are well known in theart.

According to a the proposed technology there is provided an imageprocessing device 30 for reconstructing image data from x-ray datameasured with an imaging system comprising at least one photon-countingdetector. The image processing device 30 is configured to obtain arepresentation of data measured by the photon-counting detector. Theimage processing device 30 is also configured to generate first imagedata based on a projection based first functional using a firstalgorithm, the projection based first functional being dependent on therepresentation of data. The image processing device 30 is alsoconfigured to update, based on a second functional that includes a modelof at least one physical effect not included in the projection basedfirst functional, the first image data to obtain second image data.

According to a particular embodiment of the proposed technology there isprovided an image processing device, wherein the image processing deviceis configured to generate first image data by using a first algorithm onthe first functional where the first functional is a penalized weightedleast squares functional.

According to another embodiment of the proposed technology there isprovided an image processing device, wherein the image processing deviceis configured to generate first image data by using a first algorithm onthe first functional where the first functional is a maximum likelihoodor a maximum a posteriori functional.

According to yet another embodiment of the proposed technology there isprovided an image processing device, wherein the first algorithm is aNewton method, an iterative coordinate descent method, a separablequadratic surrogates method, an expectation maximization method, aconjugate gradient method, an ordered subset method, or a combination ofa plurality of the methods.

According to still another embodiment of the proposed technology thereis provided an image processing device, wherein the processing device isconfigured to of update the first image data to obtain second image databy performing an image updating algorithm on the first functional whenat least one term from the second functional has been added to the firstfunctional, thereby yielding a perturbed first functional, the at leastone term providing a model of at least one physical effect not includedin the first functional.

By way of example, the proposed technology provides an image processingdevice that is configured to update the first image data by beingconfigured to perform a series expansion of the optimum of the perturbedfirst functional as a function of the input data, with respect to atleast one parameter determining the magnitude of at least one particularterm in the second functional which is not present in the firstfunctional, and configured to compute a correction to the first imagedata based on the series expansion, to thereby obtain second image data.

A possible embodiment of the proposed image processing device providesan image processing device wherein the second algorithm is a Newtonmethod, an iterative coordinate descent method, a separable quadraticsurrogates method, an expectation maximization method, a conjugategradient method, an ordered subset method, or a combination of aplurality of the methods.

A possible embodiment of the proposed technology provides an imageprocessing device that is configured to compute the value of the secondfunctional when the second functional:

a) depends on the first image data, and

b) depends on the second image data.

The image processing device is also configured to compare the computedvalues and select, based on the comparison, image data to be used asreconstructed image data.

Yet another embodiment of the proposed technology provides an imageprocessing device that is configured to use the selected image data asinput in an updating scheme for updating the reconstructed image data.

According to a particular example of the proposed technology there isprovided an image processing device that is configured to compute thevalue of the second functional by being configured to compute anestimate of the value of the second functional by performing at leastone step in an algorithm.

By way of example, the proposed technology provides an image processingdevice that is configured to select the image data to be used as inputby selecting the image data that yields the lowest or the highest valuefor the second functional.

A particular embodiment of the proposed technology provides an imageprocessing device wherein the physical effect is one of the following,or a combination thereof: Poisson noise statistics, optical blur,pile-up, detector element cross-talk and object scatter.

A preferred embodiment of the proposed technology provides an imageprocessing device wherein the first image data comprises at least onebasis material image.

Reference is now made to FIG. 8 where it is illustrated a schematicblock diagram of an image processing device 30, based on aprocessor-memory implementation according to an embodiment. In thisparticular example, the image processing device comprises at least oneprocessor 110 and memory 120, the memory 120 comprising instructions,which when executed by the at least one processor 110, cause the atleast one processor 110 to reconstruct image data from x-ray datameasured with an imaging system comprising at least one photon-countingdetector.

Optionally, the arrangement 100 may also include a communication circuit130. The communication circuit 130 may include functions for wiredand/or wireless communication with other devices and/or network nodes inthe network. In a particular example, the communication circuit 130 maybe based on radio circuitry for communication with one or more othernodes, including transmitting and/or receiving information. Thecommunication circuit 130 may be interconnected to the processor 110and/or memory 120. By way of example, the communication circuit 130 mayinclude any of the following: a receiver, a transmitter, a transceiver,input/output (I/O) circuitry, input port(s) and/or output port(s). FIG.8 also provides a schematic illustration of an image processing device30 that comprises a communication circuitry 130.

Alternatively, or as a complement, the arrangement may be based on ahardware circuitry implementation. Particular examples of suitablehardware circuitry include one or more suitably configured or possiblyreconfigurable electronic circuitry, e.g. Application SpecificIntegrated Circuits (ASICs), Field Programmable Gate Arrays (FPGAs), orany other hardware logic such as circuits based on discrete logic gatesand/or flip-flops interconnected to perform specialized functions inconnection with suitable registers (REG) and/or memory units (MEM).

It is also possible to provide a solution based on a combination ofhardware and software. The actual hardware-software partitioning can bedecided by a system designer based on a number of factors includingprocessing speed, cost of implementation and other requirements.

FIG. 7 is a schematic diagram illustrating an example of acomputer-implementation 200 according to an embodiment. In thisparticular example, at least some of the steps, functions, procedures,modules and/or blocks described herein are implemented in a computerprogram 225; 235, which is loaded into the memory 220 for execution byprocessing circuitry including one or more processors 210. Theprocessor(s) 210 and memory 220 are interconnected to each other toenable normal software execution. An optional input/output device 240may also be interconnected to the processor(s) 210 and/or the memory 220to enable input and/or output of relevant data such as inputparameter(s) and/or resulting output parameter(s).

The term ‘processor’ should be interpreted in a general sense as anysystem or device capable of executing program code or computer programinstructions to perform a particular processing, determining orcomputing task.

The processing circuitry including one or more processors 210 is thusconfigured to perform, when executing the computer program 225,well-defined processing tasks such as those described herein.

The processing circuitry does not have to be dedicated to only executethe above-described steps, functions, procedure and/or blocks, but mayalso execute other tasks.

In a particular embodiment, the computer program is adapted toreconstruct image data from x-ray data measured with an imaging systemcomprising at least one photon-counting detector and comprisesinstructions, which when executed by at least one processor, cause theprocessor(s) to:

-   -   read a representation of data measured by a photon-counting        detector of an x-ray imaging system;    -   generate first image data based on a first functional using a        first algorithm, the first functional being dependent on the        representation of data;    -   update, based on a second functional that includes a model of at        least one physical effect not included in the first functional,        the first image data to obtain second image data.

The proposed technology also provides a carrier comprising the computerprogram, wherein the carrier is one of an electronic signal, an opticalsignal, an electromagnetic signal, a magnetic signal, an electricsignal, a radio signal, a microwave signal, or a computer-readablestorage medium.

Hence, there is provided a computer-program product comprising acomputer-readable medium 220; 230 having stored thereon a computerprogram 225; 235.

The embodiments described above are merely given as examples, and itshould be understood that the proposed technology is not limitedthereto. It will be understood by those skilled in the art that variousmodifications, combinations and changes may be made to the embodimentswithout departing from the present scope as defined by the appendedclaims. In particular, different part solutions in the differentembodiments can be combined in other configurations, where technicallypossible.

The invention claimed is:
 1. A method of reconstructing image data fromx-ray data measured with an imaging system including at least onephoton-counting detector, the method comprising: obtaining arepresentation of data measured by said at least one photon-countingdetector, the representation comprising tomographic x-ray data measuredwith the imaging system comprising the at least one photon-countingdetector; generating first image data comprising at least one basismaterial image by optimizing a projection-based first functional using afirst optimization algorithm, said projection-based first functionalbeing dependent on said representation of data; and updating, based on asecond functional that comprises terms modelling at least one physicaleffect not included in said projection based first functional, saidfirst image data to obtain second image data, the updating comprisingperforming an image updating algorithm on the first functional when atleast one term from the second functional has been added to the firstfunctional, the at least one term providing a model of at least onephysical effect included in the second functional but not included inthe first functional, the physical effect being one or more of: opticalblur, pile-up, detector element cross-talk, and object scatter.
 2. Themethod according to claim 1, wherein the generating first image datacomprises using the first algorithm on said first functional, said firstfunctional being a penalized weighted least squares functional.
 3. Themethod according to claim 1, wherein the generating first image datacomprises using the first algorithm on said first functional, said firstfunctional being a maximum likelihood or a maximum a posteriorifunctional.
 4. The method according to claim 1, wherein said firstalgorithm is a Newton method, an iterative coordinate descent method, aseparable quadratic surrogates method, an expectation maximizationmethod, a conjugate gradient method, an ordered subset method, or acombination of a plurality of said methods.
 5. The method according toclaim 1, wherein the updating said first image data comprises performinga series expansion of the optimum of said first functional when at leastone term from said second functional has been added to said firstfunctional, said series expansion being done with respect to at leastone parameter determining the magnitude of at least one particular termfrom said second functional and computing, based on said seriesexpansion, an update of said first image data to thereby obtain secondimage data.
 6. The method according to claim 1, wherein said secondalgorithm is a Newton method, an iterative coordinate descent method, aseparable quadratic surrogates method, an expectation maximizationmethod, a conjugate gradient method, an ordered subset method, or acombination of a plurality of said methods.
 7. The method according toclaim 1, further comprising: computing the value of said secondfunctional when said second functional: a) depends on said first imagedata, and b) depends on said second image data, and comparing thecomputed values, and selecting, based on said comparison, image data tobe used as reconstructed image data.
 8. The method according to claim 7,wherein said selected image data is used as input in an updating schemefor updating said reconstructed image data.
 9. The method according toclaim 7, wherein the computing the value of said second functionalcomprises computing an estimate of the value of said second functionalby performing at least one step in an algorithm.
 10. The methodaccording to claim 8, wherein the selecting the image data to be used asinput comprises selecting the image data that yields the lowest value orthe highest value for said second functional.
 11. An image processingdevice for reconstructing image data from x-ray data measured with animaging system including at least one photon-counting detector, wherein:the image processing device is configured to obtain a representation ofdata measured by said at least one photon-counting detector, therepresentation comprising tomographic x-ray data measured with theimaging system comprising the at least one photon-counting detector; theimage processing device is configured to generate first image datacomprising at least one basis material image by optimizing aprojection-based first functional using a first optimization algorithm,said projection-based first functional being dependent on saidrepresentation of data; and the image processing device is configured toupdate, based on a second functional that comprises terms modelling atleast one physical effect not included in said projection based firstfunctional, said first image data to obtain second image data, theprocessing device being configured to update the first image data toobtain second image data by performing an image updating algorithm onthe first functional when at least one term from the second functionalhas been added to the first functional, the at least one term providinga model of at least one physical effect included in the secondfunctional but not included in the first functional, the physical effectbeing one or more of: optical blur, pile-up, detector elementcross-talk, and object scatter.
 12. The image processing deviceaccording to claim 11, wherein the image processing device is configuredto generate the first image data by using the first algorithm on saidfirst functional, said first functional being a penalized weighted leastsquares functional.
 13. The image processing device according to claim11, wherein the image processing device is configured to generate thefirst image data by using the first algorithm on said first functional,said first functional being a maximum likelihood or a maximum aposteriori functional.
 14. The image processing device according toclaim 11, wherein said first algorithm is a Newton method, an iterativecoordinate descent method, a separable quadratic surrogates method, anexpectation maximization method, a conjugate gradient method, an orderedsubset method, or a combination of a plurality of said methods.
 15. Theimage processing device according to claim 11, wherein the imageprocessing device is configured to update said first image data byperforming a series expansion of the optimum of said first functionalwhen at least one term from said second functional has been added tosaid first functional, said series expansion being done with respect toat least one parameter determining the magnitude of at least oneparticular term from said second functional, and configured to compute,based on said series expansion, an update of said first image data tothereby obtain second image data.
 16. The image processing deviceaccording to claim 11, wherein said second algorithm is a Newton method,an iterative coordinate descent method, a separable quadratic surrogatesmethod, an expectation maximization method, a conjugate gradient method,an ordered subset method, or a combination of a plurality of saidmethods.
 17. The image processing device according to claim 11, wherein:the image processing device is configured to compute the value of saidsecond functional when said second functional: a) depends on said firstimage data, and b) depends on said second image data, and the imageprocessing device is configured to compare the computed values, and theimage processing device is configured to select, based on saidcomparison, image data to be used as reconstructed image data.
 18. Theimage processing device according to claim 17, wherein the imageprocessing device is configured to use said selected image data as inputin an updating scheme for updating said reconstructed image data. 19.The image processing device method according to claim 17, wherein theimage processing device is configured to compute the value of saidsecond functional by being configured to compute an estimate of thevalue of said second functional by performing at least one step in analgorithm.
 20. The image processing device according to claim 18,wherein the image processing device is configured to select the imagedata to be used as input by selecting the image data that yields thelowest value or the highest value for said second functional.
 21. Theimage processing device according to claim 11, wherein said imageprocessing device comprises at least one processor and memory, thememory comprising instructions, which when executed by the at least oneprocessor, cause the at least one processor to reconstruct image datafrom x-ray data measured with an imaging system comprising at least onephoton-counting detector.
 22. The image processing device according toclaim 11, wherein said image processing device comprises a communicationcircuitry.
 23. A non-transitory computer program product, configured toreconstruct image data from x-ray data measured with an imaging systemincluding at least one photon-counting detector, said non-transitorycomputer program product comprising instructions, which when executed byat least one processor, cause the at least one processor to: read arepresentation of data measured by the at least one photon-countingdetector of the x-ray imaging system, the representation comprisingtomographic x-ray data measured with the imaging system comprising theat least one photon-counting detector; generate first image datacomprising at least one basis material image by optimizing a firstfunctional using a first optimization algorithm, said first functionalbeing dependent on said representation of data; and update, based on asecond functional that includes a model of at least one physical effectnot included in said first functional, said first image data to obtainsecond image data, the updating comprising performing an image updatingalgorithm on the first functional when at least one term from the secondfunctional has been added to the first functional, the at least one termproviding a model of at least one physical effect included in the secondfunctional but not included in the first functional, the physical effectbeing one or more of: optical blur, pile-up, detector elementcross-talk, and object scatter.